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Evaluation  of  a  Particle  Method  for  the  Ellipsoidal 
Statistical  Bhatnagar-Gross-Krook  Equation 

Jonathan  M.  Burt 1  and  Iain  D.  Boyd  2 

Department  of  Aerospace  Engineering 
University’  of  Michigan,  Ann  Arbor,  MI  48109 


A  particle  method  is  presented  for  the  numerical  simulation  of  rarefied  gas  flows,  based 
on  the  ellipsoidal  statistical  Bhatnagar-Gross-Krook  (ES-BGK)  model  of  the  Boltzmann 
equation.  The  simulation  procedure  includes  consideration  of  rotational  nonequilibrium, 
and  enforces  exact  momentum  and  energy  conservation  for  a  mixture  involving  monatomic 
and  diatomic  species.  This  method  is  applied  to  the  simulation  of  a  nozzle  flow  of  the  type 
associated  with  cold-gas  spacecraft  thrusters,  and  flowfield  characteristics  are  compared 
with  experimental  data  as  well  as  results  from  direct  simulation  Monte  Carlo  (DSMC)  and 
Navier-Stokes  simulations  of  the  same  flow.  The  ES-BGK  method  is  shown  to  allow  for  a 
relatively  high  degree  of  accuracy  in  transitional  flow  regimes,  while  avoiding  the 
intermolecular  collision  calculations  which  typically  make  the  DSMC  simulation  of  low 
Knudsen  number  flows  prohibitively  expensive. 


I.  Introduction 

In  the  design  and  performance  analysis  of  low-thrust  spacecraft  propulsion  systems,  various  numerical  simulation 
techniques  may  be  employed  to  determine  efficiency,  thrust  characteristics,  or  the  potential  for  plume 
impingement  and  contamination  on  spacecraft  surfaces.  A  particular  challenge  in  the  simulation  of  small  thrusters 
involving  a  chemically  inert  neutral  gas,  such  as  electro-thermal  or  cold-gas  thrusters,  is  the  accurate  consideration 
of  a  wide  range  of  Knudsen  number  regimes.  In  a  typical  thruster  of  this  type,  gas  is  expelled  through  a  convergent- 
divergent  (Laval)  nozzle  into  a  near  vacuum,  with  subsonic  near-equilibrium  flow  in  the  convergent  section  of  the 
nozzle.  As  the  gas  accelerates  through  the  divergent  section  beyond  the  nozzle  throat,  a  subsonic  boundary  layer 
grows  along  the  wall,  while  a  supersonic  core-flow  region  exists  around  the  nozzle  centerline.  The  gas  density 
continues  to  decrease  with  downstream  distance  through  the  nozzle,  and  rarefaction  effects  become  more  prominent 
within  both  the  supersonic  and  subsonic  regions.  Here  the  gas  velocity  distribution  begins  to  diverge  significantly 
from  the  equilibrium  limit,  and  thermal  energy  is  increasingly  distributed  non-uniformly  among  the  translational  and 
internal  degrees  of  freedom.  Beyond  the  nozzle  exit  plane,  these  nonequilibrium  effects  continue  to  increase  as  the 
gas  rapidly  expands  and  thermal  energy  is  converted  into  energy  associated  with  bulk  motion  of  the  exhaust  flow. 
Rotational  freezing  occurs  due  to  the  large  gradients  and  low  collision  frequency,  and  the  flow  approaches  the  free 
molecular  limit  within  a  short  distance  of  the  nozzle  exit,  particularly  at  points  far  from  the  nozzle  centerline. 

The  simulation  of  highly  rarefied  flows,  as  described  above  for  the  divergent  nozzle  region  and  plume,  is 
typically  performed  using  the  direct  simulation  Monte  Carlo  (DSMC)  method  of  Bird.1  This  method  approximates  a 
numerical  solution  to  the  Boltzmann  equation  -  the  governing  equation  for  dilute  gas  flows  based  on  a  statistical 
representation  of  molecular  velocities  -  by  decoupling  in  time  the  collision  and  advection  terms  in  the  equation.  A 
large  number  of  particles,  each  representing  a  large  number  of  atoms  or  molecules,  are  tracked  through  a 
computational  grid,  and  are  sorted  into  cells  according  to  their  location.  During  each  time  step,  some  fraction  of  the 
particles  in  a  cell  collide  with  each  other,  and  probabilistic  techniques  are  used  for  calculations  of  individual 
collisions.  All  particles  are  then  moved  through  the  grid  according  to  assigned  velocities,  and  particles  are  created  or 
removed  at  inflow  and  outflow  boundaries.  Finally,  macroscopic  quantities  are  sampled  by  averaging  various 
particle  properties  in  each  cell,  and  the  process  is  then  repeated  at  the  next  time  step. 

The  DSMC  method  has  been  shown  to  provide  accurate  solutions  for  highly  rarefied  nozzle  and  plume  flows,2 
and  retains  validity  through  the  full  range  of  possible  Knudsen  numbers  for  a  dilute  gas.  In  theory  then,  DSMC  may 
be  applied  to  the  entire  flowfield  for  a  small  spacecraft  thruster  as  described  above.  The  main  obstacle,  however,  to 
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the  DSMC  simulation  of  the  higher  density  flow  near  and  upstream  of  the  nozzle  throat  is  the  enormous 
computational  cost.  There  are  three  major  factors  which  tend  to  make  DSMC  prohibitively  expensive,  even  on  large 
parallel  machines,  for  the  simulation  of  relatively  high-density  flows:  First,  for  accurate  simulation  the  cell  size 
should  ideally  be  much  smaller  than  any  characteristic  length  scale  for  gradients  in  macroscopic  flow  quantities,  and 
nonphysical  diffusion  effects  become  increasingly  significant  as  the  cell  size  increases.  This  generally  limits  cell 
dimensions  to  roughly  one  mean  free  path,  so  that  a  very  large  number  of  cells  may  be  required  for  the  simulation  of 
even  small-scale  flows  when  the  gas  density  is  high.  Second,  the  temporal  decoupling  of  collision  and  advection 
processes  is  in  general  a  valid  approximation  only  if  the  time  step  is  no  larger  than  the  local  mean  collision  time. 
This  means  that  for  a  high  density  or  low  Knudsen  number  flow,  a  very  large  number  of  time  steps  may  be  required 
to  reach  steady  state  conditions.  These  two  requirements  may  be  relaxed  somewhat  when  macroscopic  gradient 
length  scales  are  very  large  compared  to  the  mean  free  path,  so  that  the  flow  approaches  a  near-homogeneous  state. 
However,  an  additional  limiting  factor  for  high  density  flows  cannot  be  relaxed  under  near-homogeneous  conditions: 
The  collisional  process  driving  a  flow  toward  equilibrium  is  modeled  in  DSMC  through  the  probabilistic  calculation 
of  individual  collisions  between  pairs  of  representative  particles.  The  collision  frequency  increases  with  the  gas 
density,  so  that  for  a  fixed  set  of  simulation  parameters,  a  higher  density  flow  corresponds  to  an  increase  in  the 
average  number  of  collisions  experienced  per  particle  during  each  time  step.  Under  conditions  typical  in  the 
convergent  nozzle  region  and  around  the  nozzle  throat  in  a  small  electro-thermal  or  cold-gas  spacecraft  thruster,  this 
last  requirement  can  make  DSMC  simulation  of  such  flows  prohibitively  expensive  on  all  but  the  largest  parallel 
machines. 

In  practice,  these  small-scale  exhaust  flow  simulations  may  be  performed  using  a  finite  volume  representation  of 
the  Navier-Stokes  (NS)  equations  in  the  high  density  regions  where  the  continuum  and  quasi-equilibrium 
assumptions  underlying  these  equations  are  considered  valid.  The  DSMC  method  is  then  used  only  in  the  lower 
density  regions  where  these  assumptions  break  down,  and  where  a  DSMC  flowfield  calculation  is  feasible  using 
available  computational  resources.  The  interface  between  NS  and  DSMC  regions  is  modeled  either  through  an 
uncoupled  matching  of  boundary  conditions  for  two  independent  simulations,  or  using  a  single  hybrid  code  where 
the  interface  location  may  be  varied  dynamically  during  the  simulation  through  the  evaluation  of  one  or  more 
continuum  breakdown  parameters.  For  either  approach,  the  accurate  determination  of  the  boundary  beyond  which 
the  NS  solution  is  invalid  may  be  difficult  and  expensive.  As  a  result,  there  is  generally  some  significant  tradeoff 
between  overall  simulation  accuracy  and  efficiency:  Increased  accuracy  results  from  a  more  precise  or  conservative 
determination  of  the  proper  boundary  for  the  NS  domain.  This  in  turn  requires  a  greater  computational  load,  either 
through  an  iterative  procedure  to  evaluate  the  boundary  then  update  the  overall  flowfield  solution,  or  through  an 
extension  of  the  DSMC  domain  into  regions  which  may  be  accurately  handled  by  the  NS  solver. 

This  tradeoff  between  accuracy  and  efficiency  may  to  some  extent  be  avoided  if  an  alternate  simulation  approach 
is  used  in  flowfield  regions  where  the  validity  of  the  Navier-Stokes  equations  is  questionable,  but  the  application  of 
DSMC  is  particularly  expensive.  This  alternate  method  must  be  more  efficient  than  DSMC  in  modeling  relatively 
high  density  flows,  able  to  account  for  nonequilibrium  phenomena  which  characterize  a  transition  Knudsen  number 
regime  out  of  range  of  NS  solutions,  and  in  good  quantitative  agreement  with  both  NS  and  DSMC  for  near¬ 
equilibrium  flows.  One  such  method  is  presented  here,  using  a  particle-based  probabilistic  solution  to  the  ellipsoidal 
statistical  Bhatnagar-Gross-Krook  (BGK)  model  of  the  Boltzmann  equation. 

In  recent  years,  several  authors  have  developed  particle  methods  in  which  standard  DSMC  collision  operations 
are  replaced  by  alternate  procedures  that  account  for  the  collisional  drive  toward  equilibrium  without  modeling 
individual  collisions.3"9  In  place  of  collision  calculations  on  pairs  of  representative  particles,  some  fraction  of 
particles  are  randomly  selected  during  each  time  step,  and  are  assigned  new  velocities  according  to  set  distribution 
functions.  In  several  of  these  methods,  the  collision  process  in  each  cell  within  the  computational  grid  is  based  on 
the  homogeneous  form  of  the  BGK  equation: 


d(nf)/dt  =  vn(J'-f )  (1) 

Here  f  is  the  velocity  distribution  function  for  a  simple  monatomic  gas,  /,  is  a  Maxwellian  distribution  with  the 
same  average  and  variance  as  f  n  is  the  number  density,  and  v  is  a  characteristic  relaxation  frequency.  This  is  a 
simplified  linear  model  of  the  highly  nonlinear  homogeneous  Boltzmann  equation,  where  the  detailed  physics  of  the 
collision  term  is  ignored  in  favor  of  a  Maxwell-type  molecular  model  in  which  the  entire  distribution  function  / 
decays  toward  equilibrium  at  a  velocity-independent  rate.  As  proposed  by  Macrossan6  and  independently 
implemented  by  Gallis  and  Torczynski  and  Andries  et  al.,9  the  method  may  be  modified  for  better  quantitative 
agreement  with  the  Boltzmann  equation  as  a  solution  to 
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d(nf)/dt  =  vn(fa-f) 


(2) 


where  fa  is  an  anisotropic  Gaussian  distribution  function.  This  modification  to  the  BGK  equation  was  proposed  by 
Holway10  and  Cercignani"  to  correct  for  the  nonphysical  constraint  in  the  BGK  model  of  unity  Prandtl  number,  and 
has  recently  been  shown  to  satisfy  Boltzmann’s  H-theoremfor  entropy  production.12  Following  Cercignani,  we  refer 
to  this  as  the  ellipsoidal  statistical  BGK  (ES-BGK)  model. 

In  the  first  section  of  this  paper,  we  extend  the  methods  of  Macrossan6  and  Gallis  and  Torczynskis  for  a  gas 
mixture  involving  monatomic  and  diatomic  species,  and  describe  a  procedure  to  enforce  exact  momentum  and 
energy  conservation  while  allowing  for  effects  of  rotational  nonequilibrium.  A  simple  homogeneous  test  case  is  then 
used  to  verify  the  accuracy  of  the  proposed  rotational  relaxation  model.  Next,  the  method  is  applied  to  a  nozzle  flow 
test  case  involving  a  small  cold-gas  thruster  with  N2  propellant  exhausting  into  a  vacuum,  and  results  are  compared 
with  DSMC,  NS,  and  experimental  data  for  the  same  flowfield.  We  demonstrate  relatively  good  agreement  between 
ES-BGK,  DSMC  and  experimental  results,  particularly  in  higher  density  regions  near  the  nozzle  throat,  and  show 
that  the  ES-BGK  flowfield  solution  incorporates  some  nonequilibrium  effects  which  are  also  captured  by  DSMC  but 
neglected  in  the  NS  calculation. 


II.  Numerical  Method 

As  discussed  above,  all  procedures  unrelated  to  intermolecular  collisions  in  the  ES-BGK  particle  method  are 
identical  to  those  in  DSMC.  However,  the  process  of  calculating  individual  representative  collisions  is  replaced  by  a 
procedure  which  involves  the  resampling  of  particle  properties  from  predetermined  distributions.  Consider  the  ES- 
BGK  collision  procedure  in  a  cell  containing  N  particles,  where  each  particle  may  be  assigned  to  one  of  several 
different  monatomic  or  diatomic  gas  species.  We  assume  vibrational  and  electronic  energy  modes  are  unexcited,  and 
the  variable  hard  sphere  (VHS)  collision  model1  is  used  to  approximate  the  variation  of  dynamic  viscosity  p  with 
translational  temperature  Tt.  As  a  first  step  in  the  method,  T,  is  calculated  as  a  function  of  the  average  momentum 
and  translational  energy  of  all  particles  in  the  cell: 

Tt  =i(N^l)^mC2^  -<mc)a/<m>«)  (3) 

Here  m  and  c  are  the  mass  and  velocity  respectively  of  an  individual  particle,  k  is  Boltzmann’s  constant,  the  operator 
<  >«  designates  an  average  over  all  particles  in  the  cell,  and  the  notation  x2  for  any  vector  x  is  used  to  represent  xx. 
The  coefficient  N/(N-1)  in  Eq.  (3)  corrects  for  a  statistical  error  which  appears  when  the  traditional  formula  for  Tt 
based  on  the  velocity  distribution  is  applied  to  a  discrete  set  of  particle  velocities.13 

Using  a  first  order  Chapman-Enskog  expansion  of  the  ES-BGK  equation,  it  can  be  shown  that  the  correct 
relations  for  heat  flux  and  shear  stress  are  recovered  in  a  near-equilibrium  flow  if  the  relaxation  frequency  is  given 
by 


(4) 


where  Pr  is  the  Prandtl  number  and  p  is  the  pressure.  From  the  ideal  gas  equation  of  state  and  the  VHS 
approximation  that  u  scales  with  temperature  to  some  constant  power  ox  Eq.  (4)  may  be  rewritten  as 


v  =  Pr  -  n k 


Lref 


V  Mref 


(5) 


where  Tref  is  a  reference  temperature  and  //r(.f  is  the  dynamic  viscosity  at  Tref.  For  a  simple  VHS  gas,  this  viscosity 
can  be  calculated  as  a  function  of  a>  and  the  molecular  collision  diameter  dref  at  the  reference  temperature: 


30(mkTrof  /  7r  )1//2 

4(5-2w)(7-2w)d?ef 


(6) 
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Note  that  in  most  cases  the  value  of  dref  is  itself  determined  from  experimental  measurements  of  //ref.  Eq.  (6)  is 
useful  however  in  a  hybrid  code  (discussed  below)  involving  both  DSMC  and  the  method  proposed  here,  as 
reference  collision  diameters  are  among  the  required  input  parameters  for  DSMC  collision  calculations.  A  simple 
approximate  extension  of  Eq.  (6)  to  a  gas  mixture  may  be  obtained  by  summing  the  product  of  //re[  and  the  time- 
averaged  species  mass  fraction  in  the  cell  over  all  gas  species. 

Once  the  relaxation  frequency  v  is  known,  Eq.  (2)  may  be  integrated  to  give  the  distribution  function/ after  some 
finite  time  step  At  as  a  linear  combination  of  the  initial  distribution  /(0)  and  a  corresponding  anisotropic  Gaussian 
distribution /}. 


/(At)  =  /( 0)  exp(-i/At)  +  /G  (1  -  exp(-i/At))  (7) 

From  Eq.  (7),  we  determine  the  number  of  particles  Nv  selected  for  reassignment  of  velocity  components  as 

Nv  =  floor [N(l  -  exp (—i/At))  +  R]  (8) 

where  v  is  evaluated  using  Eq.  (5),  R  is  a  random  number  in  [0,1]  and  the  operator  “floor”  rounds  to  the  nearest 
smaller  integer.  Note  that  a  time-averaged  value  of  the  number  density  n  should  be  used  in  Eq.  (5)  to  avoid  an  error 
associated  with  statistical  fluxuations  in  the  calculated  relaxation  frequency.  Once  Nv  is  known,  this  number  of 
unique  particles  is  then  randomly  selected  and  assigned  to  a  list  for  velocity  resampling.  (In  order  to  reduce  calls  to 
the  random  number  generator,  if  Nv  >  'AN  then  N  -  Nv  randomly  chosen  particles  are  removed  from  a  list  which 
initially  includes  all  particles  in  the  cell,  and  the  remaining  particles  are  used  for  resampling.)  In  the  special  case 
where  Nv=  1,  momentum  conservation  requires  that  the  velocity  of  the  selected  particle  remains  constant,  so  that  Nv 
values  of  zero  and  one  are  equivalent  during  the  resampling  procedure.  To  correct  for  the  resulting  error,  if  Nv  is 
calculated  as  one  through  Eq.  (8)  then  an  additional  random  number  is  generated  and  Nv  is  set  with  equal  probability 
to  either  zero  or  two. 

After  Nv  particles  have  been  selected,  all  velocity  components  of  these  particles  are  resampled  from  a  normalized 
Maxwellian  distribution.  Each  velocity  component  is  set  to  either 

c*  =  cos(2ttRi)^]  —ln(R2)  /  m  or  c*  =  sin  (27 tRi)^—  ln(i^)  /  m  (9) 

where  the  subscript  i  denotes  the  direction,  m  is  the  particle  mass,  and  R\  and  R2  are  random  numbers  between  0  and 
1.  Two  statistically  independent  values  of  C;  may  be  obtained  if  Eqs.  (9)  are  evaluated  for  the  same  values  of  R\  and 
R2,  so  that  only  a  single  random  number  must  be  generated  per  velocity  component.  Following  Gallis  and 
Torczynski.  we  then  modify  all  resampled  velocities  to  new  values  C;(1)  which  conform  to  an  anisotropic  Gaussian 
distribution: 


where 


/!)  _ 


=  Vi 


(10) 


Sij 


=  k- 


1  —  Pr 
2Pr 


1  N 
kTt  N  -  1 


((mCiCj 


(mci  >a(mCj  )a/<m>a)-<Sij 


(11) 


Here  </  is  the  Kronecker  delta,  c,  and  c,  are  particle  velocity  components  before  resampling,  Tt  is  evaluated  using 
Eq.  (3),  and  the  superscript  (1)  indicates  a  newly  reassigned  value. 

To  account  for  the  effect  of  collisions  on  the  rotational  energy  distribution  among  diatomic  gas  species,  we 
define  a  rotational  relaxation  frequency  vt  such  that,  in  a  homogeneous  gas,  the  instantaneous  cell-averaged 
rotational  temperature  Tr  will  vary  as 


cfll  /  dt  =  -vr  (Tr  -  Tcq )  (12) 

Here  Teq  is  the  equilibrium  temperature  to  which  both  Tr  and  Tt  converge  after  a  long  period  of  time,  and  T,  is  given 
by 
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(13) 


Tr  =^<£r)a/(Cr>a 

The  above  symbols  Ey  and  Cr  denote  rotational  energy  and  rotational  degrees  of  freedom,  respectively,  for  an 
individual  particle.  The  equilibrium  temperature  can  be  expressed  as  a  weighted  sum  of  the  rotational  and 
translational  temperatures,  such  that 


Ct,effTt  +  (Cr  )a  Tr  —  ( Ct,eff  +  (Cr  )a  )Teq  (14) 

where  Ctcfr  is  the  effective  number  of  translational  degrees  of  freedom  for  all  particles  in  the  cell.  If  all  N  particles 
are  together  considered  as  an  isolated  system,  then  it  follows  from  Eqs.  (3)  and  (13)  and  energy  conservation 
arguments  that 


d_ 

dt 


<  me  >i  /  <  m  >a  +3 


N-l\ 

— JkTt+(Cr)akTr 


=  0 


(15) 


while  momentum  and  mass  conservation  require  that  the  first  term  within  the  brackets  be  constant  in  time  as  well. 
By  comparing  Eqs.  (14)  and  (15)  and  noting  that  Teq  is  a  time  invariant  quantity,  we  can  show  that  there  are 
effectively  3(N-1)/N  translational  degrees  of  freedom  in  the  cell,  so  that  the  equilibrium  temperature  may  be  given 
from  Eq.  (14)  as 


_  3((N-l)/N)Tt  +(Cr)aTr 

eq  3(N-1)/N  +  (Cr)a  1  J 

Note  that  the  arguments  used  to  derive  Eq.  (16)  make  no  assumptions  regarding  how  the  translational  and 
rotational  temperatures  relax  toward  the  equilibrium  state,  so  that  this  equation  should  also  be  valid  for  DSMC  and 
other  methods  for  which  a  discrete  set  of  particle  velocities  is  used  to  describe  the  velocity  distribution  of  a  gas  in 
rotational  nonequilibrium.  As  there  are  exactly  three  translational  degrees  of  freedom  for  any  real  gas  under 
consideration  here,  this  derivation  reveals  the  existence  of  a  deterministic  error  which  scales  as  (N-l)/N  and  may 
significantly  influence  calculated  bulk  flow  properties  when  the  average  number  of  particles  per  cell  is  very  small. 
To  the  authors’  knowledge  no  previous  references  to  this  error  source  exist  in  the  literature,  although  the  upper 
bound  of  the  resulting  error  should  be  relatively  small  (less  than  5%)  when  particle  numerical  weights  are  set  such 
that  over  20  particles  are  assigned  to  each  cell,  as  is  standard  practice  in  DSMC. 

The  rotational  relaxation  frequency  vT  shown  in  Eq.  (12)  is  defined  here  as  the  frequency  of  inelastic  collisions, 
given  as  the  ratio  of  the  collision  frequency  Fcoll  to  the  rotational  collision  number  Zr.  To  determine  Fcon,  consider  a 
cell  containing  particles  of  Nspec  different  species,  each  with  a  mass  m,.  VHS  reference  collision  diameter  drcr,  and 
number  density  /?,  averaged  over  a  large  number  of  time  steps.  Assuming  a  near-equilibrium  flow,  we  can  efficiently 
calculate  the  collision  frequency  as 


^coll 

n 


/kTref 

f  Tt  ) 

1  27 r 

,  Tref , 

1 — uj  Nspec  i 

>2  !EninjAij(2  -  h) 

i=i  j=i 


where 

A  TT/j  ,  \2  lm  i  +  mi  , 

Aij  —  4  (dref.i  +  ^ref  j  J  J  an^  11  ~  ni  • 

V  i  AAj  i=l 

Note  that  no  implied  summations  are  used  in  the  index  notation  of  Eq.  (17). 

The  rotational  collision  number  Zr  is  calculated  using  a  modified  form  of  a  formula  given  by  Parker.14 


_ 3/5  Zg° _ 

1  +  (tt1/2  /2)(T*/Teq)1'2  +  (tt  +  Tr2  /4)(T*/Teq) 


(17) 


(18) 
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The  above  constants  Zrx  and  T  are  determined  for  N2  as  Zrx=23.0  and  T*=91.5  K.1  As  a  modification  to  Parker’s 
formula,  the  equilibrium  temperature  Teq  is  used  here  in  place  of  Tt  following  observations  of  Boyd15  that,  for  a 
single  intermolecular  collision,  Z,  tends  to  vary  as  a  function  of  both  the  translational  and  rotational  collision 
energies.  The  3/5  factor  in  the  numerator,  also  not  in  Parker’s  original  formula,  accounts  for  the  fact  that  Z,  is 
defined  here  to  characterize  the  rate  at  which  Tr  relaxes  toward  Teq,  not  toward  Tt  as  used  by  Parker  and  others. 1 

By  analogy  with  Eq.  (8),  Eq.  (12)  can  be  integrated  in  time  to  give  an  initial  estimate  of  the  number  of  particles 
N,  for  which  to  resample  rotational  energy  as 


Nr 


N 


1  —  exp 


(19) 


where  Fcan  and  Z,  are  evaluated  using  Eqs.  (17)  and  (18)  respectively.  (The  exact  determination  of  Nr  requires  a 
more  involved  procedure,  and  is  discussed  in  detail  below.)  In  our  implementation,  Nr  particles  are  randomly 
selected  from  the  list  of  Nv  particles  which  have  already  been  tagged  for  velocity  resampling,  and  the  rotational 
energy  Et  of  each  selected  diatomic  particle  is  reassigned  to  a  normalized  value 

4 D  =  —  ln(i?)  (20) 


where  a  unique  random  number  R  e  [0,1]  is  used  for  each  particle. 

To  avoid  random  walk  errors  and  permit  overall  accuracy  for  a  near-equilibrium  flow  comparable  to  that  of 
either  DSMC  or  NS  simulations,  we  desire  for  both  momentum  and  energy  to  be  exactly  conserved  during  the 
resampling  procedures.  It  follows  that  all  reassigned  velocities  c(1)  and  rotational  energies  1 1  must  be  modified  to 
final  values  c<2)  and  Er{1]  such  that 


|mc,]|  |  =  |mc,0'|  (21) 

and 

iNv  (m(c(2))2  )y  +  Nr  (42)  )r  =  ^Nv  (m(c(°))2  )y  +  Nr  (40)  )r  (22) 

where  the  operators  <  >v  and  <  >r  denote  averages  over  all  particles  used  in  the  resampling  of  velocity  and 
rotational  energy  respectively.  Further,  to  satisfy  energy  equipartition  at  the  equilibrium  limit,  we  require  that 
translational  and  rotational  temperatures  based  on  reassigned  values  be  equal.  This  condition  can  be  written  as 

)y  -  (mc(2)  )v  /<m>V  )  =  ^(^2))r/<Cr)r  (23) 

where  C,T  is  the  number  of  rotational  degrees  of  freedom  for  an  individual  particle.  For  each  of  the  Nv  particles 
selected  for  velocity  resampling,  the  final  velocity  is  calculated  as 

=  B  +  C  (c(1)  —  ^mc*1)  ^  /<m>v  )  (24) 

while  each  of  the  Nr  particles  with  resampled  rotational  energies  is  given  a  final  value  of 

Ef>=DE 4.  (25) 


By  algebraic  manipulation  of  Eqs.  (21)  through  (25),  we  can  show  that  momentum  conservation,  energy 
conservation  and  energy  equipartition  are  all  satisfied  if  the  vector  B  and  scalar  constants  C  and  D  are  chosen  as 
follows: 


B  =  ^  me®  ^  /  <  m  >v 


(26) 
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,1/2 


(27) 


f  3(Nv  —  1)  ] 

^m(c1®)2  ^  —  ^mc(0)  ^  /<m)v  +  2(Nr  /Nv  ){E^  ^ 

l3(Nv-l)  +  Nr(C)rJ 

^m(c(1))2^  -|mc^|2/(m)v 

(  Nv(C)r  ] 

^  m(c(0)  )2  ^  ^  me®  ^  /  <  m  >v  +  2  ( Nr  /  Nv )  ^  ^ 

(3(Nv-l)  +  Nr<C)rJ 

2(Et«)r 

(28) 


From  the  definition  of  rotational  temperature  indicated  by  Eq.  (13),  the  right  hand  side  of  Eq.  (28)  may  be 
expressed  as  a  ratio  of  temperatures:  The  denominator  gives  the  rotational  temperature  based  on  values  of  E:(  1 \ 
while  the  numerator  gives  the  temperature  corresponding  to  the  conditions  of  Eq.  (23).  Thus,  Eq.  (28)  may  be  used 
to  show  that  the  temperature  T<2)  based  on  resampled  and  renormalized  values  of  either  particle  velocity  c<21  or 
rotational  energy  Ep]  can  be  given  as 


t(2) 


V 


3(NV  —  l)  +  Nr(Cr)r 


[mc,0))  /<m>v 


vNv 


(29) 


As  particles  are  randomly  chosen  for  velocity  or  rotational  energy  resampling,  the  averaged  quantities  in  Eq.  (29) 
will,  on  average,  equal  the  corresponding  quantities  averaged  over  all  N  particles  in  the  cell.  If  we  denote  the 
equivalent  temperature  based  on  cell-averaged  values  as  Ta<2),  then  it  follows  from  substitution  of  Eqs.  (3)  and  (13) 
into  Eq.  (29)  that 


t(2)  =  3((Nv-l)/Nr)Tt+(Cr)aTr 
a  3(Nv-l)/Nr+(C)a  1  J 

As  the  temperature  Ta<2)  is  generally  not  equal  to  the  equilibrium  temperature  Teq  determined  from  Eq.  (16),  the 
procedure  used  to  calculate  Nr  must  be  based  on  a  relaxation  frequency  other  than  the  true  characteristic  frequency 
of  rotational  relaxation  vr  =  Fc 0u/  Zr.  A  correction  factor  J3  is  therefore  introduced  as  the  ratio  of  this  required 
relaxation  frequency  to  vr,  so  that  in  a  homogeneous  adiabatic  flow, 


c/T  /  dt  =  -pvr  (T  -  Ta,2) ) 


From  a  comparison  of  Eqs.  (12)  and  (31),  it  follows  that 

A  =  (Tr-Teq)/(Tr-Tf) 


Substituting  Eqs.  (16)  and  (30)  into  Eq.  (32),  we  find 


3  +  (  Cr  )a  Nr*/(NV-1) 
3  +  (Cr  )a  N/(N-1) 


(31) 


(32) 


(33) 


where  N/  is  an  approximation  to  the  integer  value  Nr  based  on  an  exact  solution  to  Eq.  (31).  If  Eq.  (31)  is  integrated 
over  the  time  step  interval  At,  then  the  same  reasoning  used  to  find  Nv  by  Eq.  (8)  may  be  used  to  calculate  N,  as 


Nr  =  N 


1  —  exp 


Zr 


(34) 


Note  from  Eqs.  (33)  and  (34)  that  [i  and  Nr*  cannot  be  independently  determined,  but  must  instead  be  calculated 
simultaneously  using  an  iterative  method.  In  the  implementation  proposed  here,  P  is  initially  set  to  one  and  N,  is 
found  by  Eq.  (34),  after  which  P  and  Nr  are  updated  using  Eqs.  (33)  and  (34)  respectively.  This  procedure  is 
repeated  four  additional  times,  after  which  the  number  of  particles  to  select  for  rotational  energy  resampling  is 
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computed  as  Nr  =  floor[N,*  +/?]  where  R  is  a  random  number  in  [0, 1].  This  simple  “brnte  force”  approach  is  found  to 
give  a  reasonable  balance  of  accuracy  and  efficiency  over  a  wide  range  of  conditions,  and  more  complicated 
iterative  schemes  with  faster  convergence  should  not  generally  be  required  due  to  the  insignificant  contribution  of 
this  procedure  to  the  total  computational  expense. 

III.  Homogeneous  Rotational  Relaxation 

To  verify  the  accuracy  of  the  procedure  described  above  for  rotational  relaxation,  this  procedure  is  implemented 
in  a  modified  version  of  the  DSMC  code  MONACO16  and  applied  to  a  zero-dimensional  test  case  based  on 
simulations  of  Bird.1  An  isolated  control  volume  of  N2  gas  is  given  initial  conditions  involving  translational 
equilibrium  at  500  K  and  a  rotational  temperature  of  0  K.  The  simulation  domain  consists  of  a  single  cell  with 
specularly  reflecting  walls  and  10,000  particles,  and  the  time  step  is  set  to  1/2  of  the  initial  mean  collision  time.  The 
rotational  collision  number  is  found  from  Eqs.  (16)  and  (18)  to  hold  a  constant  value  of  Zr  =  4.31,  so  that  the 
rotational  and  translational  temperatures  will  vary  in  time  according  to  analytical  expressions  similar  to  those 
derived  by  Bird.1 


Tr  =  300 


1  —  exp 


4.31  rcon 


,  Tt  =  300  +  200  exp 


4.31  rcon 


(35) 


Here  both  temperatures  are  given  in  K,  t  is  the  total  elapsed  time  and  Tcon  is  a  characteristic  mean  collision  time.  The 
variable  xcon  may  be  defined  in  terms  of  the  instantaneous  mean  collision  frequency  Fcon  as 


Tcoll 


(36) 


where  Fa,n  is  calculated  for  a  VHS  gas  using  Eq.  (17). 

Figure  (1)  shows  the  variation  in  both  Tr  and  Tt  as  functions  of  the  normalized  time  t/icon.  Results  from  the  ES- 
BGK  simulation  are  presented  along  with  results  from  a  DSMC  simulation  under  identical  conditions,  and 
theoretical  curves  based  on  Eqs.  (35)  are  also  shown  for  comparison.  For  consistency,  a  constant  rotational  collision 
number  of  4.31  is  used  in  the  DSMC  simulation.  The  ES-BGK  data  points  are  found  to  be  in  close  agreement  with 
both  theory  and  DSMC,  while  a  slight  underprediction  is  observed  in  the  relaxation  rate  from  the  ES-BGK 
simulation.  This  discrepancy  can  be  traced  to  the  fact  that  only  a  fraction  of  particles  are  chosen  for  velocity 
resampling  during  each  time  step.  If  the  artificial  condition  NV=N  is  imposed,  then  nearly  perfect  agreement  is  found 
between  the  ES-BGK  and  theoretical  curves  of  Fig.  (1).  As  discussed  by  Bird,1  even  DSMC  results  for 
homogeneous  rotational  relaxation  with  constant  Z,  can  only  be  approximated  by  an  analytical  solution,  and  the 
scale  of  the  disagreement  between  DSMC  and  theoretical  curves  is  comparable  to  that  between  theory  and  ES-BGK. 
Both  DSMC  and  ES-BGK  results  have  some  dependence  on  the  VHS  temperature  exponent  to,  so  by  comparison 
with  DSMC  we  can  conclude  that  rotational  relaxation  is  modeled  with  acceptable  accuracy  by  the  method  proposed 
here. 


IV.  Nozzle  Flow  Simulation 

Few  experimental  studies  are  available  in  the  literature  involving  the  measurement  of  flowfield  characteristics 
within  a  nozzle  of  the  type  employed  in  cold-gas  spacecraft  thrusters.  To  the  authors’  knowledge,  the  most  detailed 
is  that  of  Rothe,17  where  an  electron  beam  technique  was  used  to  measure  density  and  rotational  temperature  for  the 
flow  of  N2  through  a  very  small  nozzle  into  a  vacuum  chamber.  Both  NS  and  DSMC  simulations  of  the  internal 
nozzle  flow  have  been  performed  by  Chung  et  al.18  In  the  DSMC  simulation  by  these  authors,  the  computational 
domain  is  restricted  to  the  lower  density  regions  beyond  the  nozzle  throat,  with  inflow  conditions  at  the  throat  based 
on  a  NS  simulation  of  the  flow  within  the  convergent  section  of  the  nozzle.  Generally  good  agreement  is  found 
between  DSMC  results  and  experimental  data,  while  a  NS  simulation  of  the  entire  nozzle  flow  gives  relatively  poor 
agreement  with  the  experiment  in  areas  downstream  of  the  throat  where  nonequilibrium  effects  become  significant. 

To  evaluate  the  accuracy  of  the  ES-BGK  approach  presented  here,  we  perform  ES-BGK  and  DSMC  simulations 
for  the  divergent  nozzle  and  nearfield  plume  flow  investigated  by  Rothe.17  Identical  grids,  boundary  conditions  and 
numerical  weights  are  used  for  the  two  simulations,  with  the  same  grid  geometry  and  inflow  conditions  at  the  nozzle 
throat  employed  in  the  DSMC  simulation  of  Chung  et  al.18  As  the  ES-BGK  particle  method  is  based  on  assumptions 
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of  near-equilibrium  flow,  we  can  investigate  the  applicability  of  the  method  to  higher  Knudsen  number  regimes  by 
comparing  ES-BGK  results  with  both  DSMC  and  experimental  data. 

Following  Chung  et  al.  we  use  a  VHS  temperature  exponent  ( oj)  of  0.24  and  assume  fully  diffuse  reflection  for 
particle  collisions  with  the  nozzle  wall.  These  authors  investigated  a  number  of  different  models  for  translational 
energy  transfer  during  particle-wall  collisions,  and  found  best  agreement  with  experimental  data  when  all  collisions 
are  assumed  diffuse,  with  a  10%  probability  of  thermal  accommodation  to  the  wall  temperature  of  300  K.  Energy 
exchange  between  translational  and  rotational  modes  was  calculated  using  a  detailed  model  which  allows  the 
exchange  probability  to  vary  as  a  function  of  the  total  incident  energy  of  the  colliding  particle.  As  noted  by  these 
authors,  a  lack  of  experimental  measurements  on  the  interaction  of  N2  gas  with  solid  graphite  (the  nozzle  wall 
material)  leaves  some  doubt  as  to  the  validity  of  this  energy  exchange  model.  A  more  commonly  accepted  procedure 
for  modeling  rotational  energy  transfer  during  particle-wall  collisions  in  DSMC  involves  an  assumption  that  particle 
rotational  energy  is  resampled  from  an  equilibrium  distribution  at  the  wall  temperature  with  some  uniform 
probability.  Therefore,  as  a  simplification  to  the  procedure  described  by  Chung  et  al.,  we  assume  that  10%  of 
particle-wall  collisions  involve  full  accommodation  of  both  translational  and  rotational  energy  to  the  wall 
temperature,  and  the  remaining  90%  of  these  collisions  involve  adiabatic  diffuse  reflection.  In  the  latter  type  of 
collision,  the  post-collision  direction  of  the  particle  trajectory  is  randomly  assigned  according  to  the  equilibrium 
velocity  distribution  for  particles  which  pass  through  a  plane,  while  the  incident  values  of  both  translational  and 
rotational  energy  are  retained. 

Axisymmetric  ES-BGK  and  DSMC  simulations  are  performed  on  a  1.4  GHz  AMD  Athlon  cluster,  with  static 
domain  decomposition  among  eight  separate  tasks.  Each  calculation  is  run  for  24  hours,  with  roughly  15  million 
particles  at  steady  state  and  an  unstructured  grid  of  about  41,000  triangular  cells  scaled  according  to  the  local  mean 
free  path.  The  grid  geometry  is  shown  in  Fig.  (2)  following  that  given  by  Chung  et  al.  Not  visible  in  the  figure  is  the 
wall  curvature  at  the  throat,  where  the  radius  of  curvature  is  one  half  the  throat  radius  of  2.55  mm.  The  nozzle 
divergence  angle  is  20°,  the  exit-to-throat  area  ratio  is  66,  and  the  simulation  domain  extends  14  mm  downstream  of 
the  exit  plane.  To  avoid  unphysical  flow  characteristics  beyond  the  exit  plane  due  to  the  small  domain  size,  we 
neglect  the  effect  of  ambient  pressure  in  the  vacuum  chamber  and  use  outflow  boundaries  within  the  plume  region. 
Results  for  both  ES-BGK  and  DSMC  simulations  are  shown  in  Figs.  (3)  through  (10). 

V.  Simulation  Results 

Mach  number  contours  for  both  simulations  are  shown  in  Fig.  (3),  where  the  top  half  and  bottom  half  of  the 
figure  give  ES-BGK  and  DSMC  results,  respectively.  Excellent  agreement  between  the  two  methods  is  found  in  the 
area  around  the  nozzle  throat.  In  both  simulations,  a  significant  reduction  in  Mach  number  is  found  toward  the 
centerline  in  this  region,  as  results  from  a  weak  compression  wave  due  to  the  curvature  of  streamlines  at  the  throat.18 
Mach  number  contours  along  the  nozzle  wall  are  also  nearly  identical,  and  both  results  show  a  termination  of  the 
sonic  line  at  the  nozzle  lip  due  to  the  rapid  acceleration  expected  at  the  lip  in  any  highly  underexpanded  nozzle  flow. 
However,  prominent  differences  between  the  two  results  appear  along  the  centerline  beyond  a  point  about  1  cm 
downstream  of  the  throat.  As  the  degree  of  rarefaction  increases  with  downstream  distance,  the  centerline  Mach 
number  gradient  in  the  DSMC  simulation  becomes  significantly  greater  than  in  the  ES-BGK  simulation,  so  that  the 
intersection  of  the  Mach  4  contour  line  with  the  axis  occurs  roughly  2.6  cm  further  upstream  in  DSMC  than  in  ES- 
BGK  results.  This  difference  in  Mach  number  is  likely  a  result  of  the  unphysical  lack  of  velocity  dependence  in  the 
collisional  relaxation  rate  of  Eq.  (2),  which  should  correspond  to  some  reduction  in  the  streamwise  temperature 
gradient  within  areas  of  low  collision  frequency.  The  difference  between  the  results  is  however  largely  restricted  to  a 
narrow  region  along  the  axis,  and  very  good  agreement  between  the  two  sets  of  contours  is  still  found  well 
downstream  of  the  throat  at  points  further  from  the  centerline.  Note  that  no  contour  lines  are  shown  within  the  plume 
backflow  region  beyond  the  nozzle  lip,  due  to  the  small  sample  size  in  this  very  low  density  region. 

Contours  of  density  are  shown  in  Fig.  (4),  again  with  ES-BGK  results  in  the  top  half  and  DSMC  results  in  the 
bottom  half  of  the  figure.  All  values  are  normalized  by  the  stagnation  density  p0  =  5.323xl0"3  kg/m'.  As  in  Fig.  (3) 
we  find  excellent  agreement  around  the  nozzle  throat  and  further  downstream  in  regions  far  from  the  centerline. 
Both  results  show  an  off-axis  local  maximum  in  the  density  profile  just  downstream  of  the  throat,  which  is 
consistent  with  a  weak  compression  wave  described  by  Chung  et  al.18  and  mentioned  above  in  the  discussion  of 
Mach  number  contours.  As  above,  the  greatest  difference  between  the  two  results  is  found  along  the  axis  at  points 
far  from  the  throat.  Here  there  is  a  noticeable  increase  in  the  ES-BGK  expansion  rate  relative  to  that  of  DSMC,  so 
that  at  a  given  point  on  the  axis  the  density  is  lower  in  the  ES-BGK  results. 

This  trend  is  also  found  in  Fig.  (5),  which  shows  the  normalized  density  variation  along  the  centerline.  The  axial 
coordinate  x  in  the  plot  is  normalized  by  the  throat  radius  Rt  =  2.55  mm,  and  a  dotted  vertical  line  gives  the  location 
of  the  nozzle  exit  plane.  Here  results  for  the  ES-BGK  and  DSMC  simulations  are  plotted  against  those  for  DSMC 
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and  NS  simulations  of  Chung  et  al. 1  and  the  experimental  data  points  of  Rothe.17  For  x/Rt  <  10  very  similar 
centerline  density  profiles  are  found  for  the  two  simulations  presented  here,  and  a  slight  difference  is  observed 
between  these  and  the  DSMC  and  NS  contours  given  by  Chung  et  al.  The  reason  for  this  discrepancy  is  not  clear, 
but  two  likely  possibilities  are  some  difference  in  grid  resolution  within  this  region  between  the  two  sets  of 
simulations,  or  a  discretization  of  flow  conditions  at  the  nozzle  throat  from  the  use  of  a  finite  number  of  uniform 
inflow  boundaries.  In  general,  the  ES-BGK  result  agrees  very  well  with  both  DSMC  density  profiles  and  the 
experimental  values.  We  find  a  significantly  higher  overall  accuracy  here  for  ES-BGK  than  for  NS,  particularly  far 
downstream  of  the  throat  where  density  is  slightly  underpredicted  by  ES-BGK  but  greatly  overpredicted  in  the  NS 
simulation. 

Figure  (6)  shows  the  variation  in  translational  and  rotational  temperature  along  the  nozzle  centerline  for  ES-BGK 
and  both  DSMC  simulations,  along  with  the  NS  centerline  temperature  and  rotational  temperature  values  of  Rothe. 
All  temperatures  here  are  normalized  by  the  stagnation  temperature  T0  =  300  K.  In  comparing  the  two  temperature 
profiles  from  the  ES-BGK  simulation,  note  that  the  rotational  temperature  is  uniformly  higher  than  the  translational 
temperature,  and  that  the  difference  between  the  two  values  increases  with  downstream  distance.  This  trend  is  also 
observed  in  both  sets  of  DSMC  results,  and  demonstrates  the  rotational  energy  lag  expected  in  rarefied  expansion 
flows  due  to  a  finite  rotational  relaxation  time.  The  quantitative  level  of  agreement  between  the  experimental 
measurements  and  ES-BGK  results  is  in  general  comparable  to  that  between  the  experiment  and  DSMC,  although 
the  ES-BGK  rotational  temperature  near  the  nozzle  exit  is  significantly  higher  than  that  found  using  the  other 
methods. 

Density  profiles  across  a  radial  plane  near  the  nozzle  exit  are  shown  in  Fig.  (7).  This  plane  is  defined  by  the  axial 
location  x/Rt  =  18.7,  and  radial  coordinates  y  are  normalized  by  the  local  nozzle  radius  Rw.  As  discussed  above,  the 
correlation  between  ES-BGK  and  DSMC  results  in  this  region  decreases  significantly  toward  the  nozzle  centerline, 
where  the  ES-BGK  density  is  about  7%  lower  than  that  found  for  either  DSMC  simulation.  Further  from  the 
centerline,  the  agreement  between  ES-BGK  and  DSMC  results  increases  considerably,  with  nearly  identical  values 
found  at  points  close  to  the  nozzle  wall.  At  all  points  along  this  plane,  ES-BGK  density  values  are  uniformly  lower 
than  those  found  experimentally,  while  the  density  values  extracted  from  the  NS  solution  are  significantly  greater 
than  the  experimental  data.  Except  for  a  small  region  roughly  half  way  between  the  axis  and  nozzle  wall,  we  find  the 
ES-BGK  density  values  to  more  closely  approximate  measured  values  than  do  results  from  the  NS  simulation. 

Profiles  of  rotational  and  translational  temperature  along  the  same  radial  plane  are  shown  in  Fig.  (8).  In 
comparing  the  ES-BGK  rotational  temperature  profile  with  that  from  the  DSMC  simulation  performed  here,  we  find 
an  overestimate  of  about  13%  in  the  ES-BGK  values  near  the  axis,  and  a  slight  underestimate  (less  than  1%)  at  the 
nozzle  wall.  The  same  general  pattern  is  also  observed  in  a  comparison  of  translational  temperature  profiles, 
although  the  translational  temperature  is  underpredicted  in  the  ES-BGK  simulation  by  roughly  3%  at  the  wall.  While 
some  of  these  discrepancies  are  significant,  the  general  trends  shown  in  the  DSMC  results  are  qualitatively 
reproduced  very  well  by  the  ES-BGK  method:  The  rotational  temperature  is  slightly  higher  than  the  translational 
temperature  for  y/Rw<  0.7,  with  a  sign  reversal  in  the  temperature  difference  beyond  this  point.  At  the  nozzle  wall, 
the  translational  temperature  is  almost  fully  accommodated  to  the  wall  temperature,  while  the  rotational  temperature 
is  considerably  lower. 

The  difference  between  translational  and  rotational  temperatures  at  the  wall  in  both  the  ES-BGK  and  DSMC 
results  may  be  attributed  to  the  wall  collision  model.  All  particle-wall  collisions  are  diffuse,  so  that  the  translational 
component  of  average  thermal  energy  for  reflected  particles  corresponds  to  either  the  wall  temperature  or  the  local 
stagnation  temperature  as  measured  by  the  velocity  of  incident  particles.  As  the  wall  and  stagnation  temperatures  are 
nearly  equal,  we  expect  a  very  small  difference  between  the  wall  surface  temperature  and  the  gas  translational 
temperature  at  the  wall.  In  contrast,  the  average  post-collision  rotational  energy  of  particles  involved  in  adiabatic 
diffuse  reflection  corresponds  to  the  rotational  component  of  the  gas  static  temperature,  not  the  stagnation 
temperature.  Rotational  energy  accommodation  occurs  for  only  a  small  fraction  of  wall  collisions,  so  there  should  be 
relatively  little  difference  in  the  average  rotational  energy  between  incident  and  reflected  particles.  In  low  density 
regions  where  the  rotational  relaxation  time  is  large  compared  to  other  characteristic  time  scales  for  the  flow,  a 
significant  lag  may  develop  between  rotational  and  translational  temperatures.  As  the  region  under  consideration  has 
both  a  large  rotational  relaxation  time  and  a  boundary  condition  which  tends  to  disproportionately  increase  the 
translational  temperature,  the  trends  observed  at  the  wall  in  Fig.  (8)  seem  physically  justified. 

Note  that  the  NS  simulation  does  not  account  for  rotational  nonequilibrium,  so  only  a  single  temperature  for  this 
simulation  is  shown  in  the  figure.  The  NS  temperature  profile  shows  relatively  poor  agreement  with  all  of  the  other 
results,  and  gives  a  gas  temperature  at  the  wall  between  the  translational  and  rotational  temperatures  found  using  the 
other  methods.  A  no-slip  adiabatic  boundary  condition  is  used  for  the  wall  in  the  NS  calculation,  so  some 
temperature  discontinuity  at  the  wall  is  expected. 
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Profiles  of  bulk  axial  velocity  are  shown  in  Fig.  (9),  again  at  the  radial  plane  near  the  nozzle  exit  defined  by  x/Rt 
=  18.7.  Velocities  here  are  normalized  by  the  most  probable  thermal  speed  at  the  stagnation  temperature  SG  = 
(2RT0)'/2  =  422  m/s.  We  find  nearly  identical  profiles  for  the  ES-BGK  and  the  two  DSMC  simulations,  with  a 
deviation  between  the  three  results  of  less  than  5%  at  all  points  on  the  plot.  The  NS  velocity  profile  shows  excellent 
agreement  with  the  other  results  at  the  centerline,  but  diverges  significantly  within  a  short  distance  from  the  axis. 
While  the  ES-BGK  and  DSMC  simulations  give  a  wall  slip  velocity  equal  to  about  28%  of  the  centerline  velocity, 
the  NS  no-slip  boundary  condition  results  in  a  velocity  of  zero  at  the  wall.  If  a  Maxwell-type  slip  model  (based  on 
the  velocity  gradient  at  the  wall)  had  been  employed  in  the  NS  simulation,  we  expect  that  better  agreement  would  be 
found  between  the  NS  velocity  profile  and  those  for  the  other  methods.  A  slip  model  would  however  affect  density 
and  temperature  profiles  as  well,  so  it  is  not  clear  that  its  inclusion  would  lead  to  a  net  gain  in  overall  accuracy  for 
the  NS  simulation  of  this  flow. 

Figure  (10)  displays  the  variation  along  the  centerline  of  continuum  breakdown  parameters,  based  on  the  DSMC 
simulation  performed  for  this  study.  As  defined  by  Boyd  et  al., 19  the  local  degree  of  thermal  nonequilibrium  may  be 
roughly  characterized  by  the  maximum  value  of  a  “gradient  length  local”  Knudsen  number  KnGLL(Q),  where 


KnGLL(Q) 


_A  dQ 

Q  ds 


(37) 


and  Q  is  set  to  either  the  density,  translational  temperature  or  bulk  speed  of  the  gas.  The  symbol  X  here  denotes  the 
mean  free  path,  and  cQ/os  is  the  partial  derivative  of  Q  in  the  direction  of  maximum  gradient.  As  radial  derivatives 
of  macroscopic  flow  quantities  are  by  definition  zero  along  the  centerline,  0s  here  is  in  the  axial  direction. 
Significant  scatter  is  observed  among  data  points  in  the  figure,  as  a  result  of  both  statistical  noise  and  subtractive 
cancellation  error  in  the  evaluation  of  the  derivatives.  Tenth  order  polynomial  trendlines,  based  on  a  least-squares 
fit,  are  therefore  used  to  more  clearly  show  the  general  trends  among  the  three  different  parameters. 

As  is  generally  expected  for  a  near-adiabatic  expansion  flow,  the  parameter  based  on  density  KnGLL(/>)  gives  the 
largest  value  through  nearly  the  entire  range  considered.  Based  on  this  parameter,  we  find  a  local  maximum  in  the 
degree  of  nonequilibrium  a  short  distance  downstream  of  the  throat.  This  maximum  can  be  attributed  to  the  large 
curvature  of  streamlines  around  the  throat,  and  the  associated  weak  compression  wave  mentioned  in  the  discussion 
of  Figs.  (3)  and  (4).  The  flow  then  equilibrates  somewhat  before  rarefaction  effects  again  become  more  prominent 
further  downstream.  A  maximum  KnGLL(/>)  value  of  slightly  less  than  0.025  is  ultimately  reached  in  the  nearfield 
plume  region  just  upstream  of  the  outflow  boundary.  Trends  in  this  parameter  are  found  to  correlate  well  with  the 
general  level  of  agreement  between  ES-BGK  and  DSMC  results  as  shown  in  Figs.  (3)  through  (6).  As  the  value  of 
KnGLL(/>)  increases  with  downstream  distance  along  the  centerline,  translational  and  rotational  temperatures  are 
increasingly  overpredicted  by  the  ES-BGK  method,  while  density  and  Mach  number  are  underpredicted. 

Note  that  the  comparison  of  ES-BGK  simulation  results  with  DSMC,  NS  and  experimental  data  is  intended  to 
assess  the  validity  of  the  assumptions  and  approximations  behind  the  method  itself,  and  not  to  evaluate  the  tradeoff 
between  efficiency  and  accuracy  at  larger  time  steps  where  advantages  of  the  ES-BGK  method  over  DSMC  become 
apparent.  Further  study  of  time  discretization  errors  on  a  simpler  transition-regime  flow  is  warranted  to  more  clearly 
show  the  potential  advantages  of  this  method.  However,  if  computational  expense  is  measured  as  the  mean 
calculation  time  per  particle  per  time  step  at  steady  state,  the  ES-BGK  simulation  presented  here  is  about  9%  less 
expensive  than  the  DSMC  simulation.  While  this  figure  may  seem  discouraging,  no  optimization  for  efficiency  is 
attempted  here  through  either  local  or  global  variation  in  the  time  step  size.  Instead,  a  single  time  step  size  is  used 
throughout  the  grid  domain  following  the  minimum  mean  collision  time  at  the  nozzle  throat.  As  a  result  of  the  large 
ratio  of  mean  collision  time  to  time  step  size  through  much  of  the  flowfield,  a  relatively  small  fraction  of  the  total 
simulation  time  is  taken  up  by  collision  operations.  Given  this  fact,  a  9%  increase  in  overall  efficiency  does  make 
the  proposed  method  seem  a  promising  alternative  to  DSMC  for  flows  involving  moderate  rarefaction  effects  where 
collision  operations  account  for  a  larger  fraction  of  the  total  simulation  time. 

As  the  time  step  and  cell  size  increase,  some  additional  loss  of  accuracy  is  expected  in  the  ES-BGK  results,  due 
to  an  error  associated  with  the  temporal  and  spatial  decoupling  of  collision  and  advection  processes  in  the  gas,  as 
well  as  a  nonphysical  equilibration  of  flow  properties  in  each  cell  and  an  increase  in  numerical  diffusion  effects.20  At 
the  same  time,  the  computational  expense  of  collision  procedures  will  approach  an  asymptotic  limit  where  all 
particle  velocity  and  rotational  energy  values  are  reassigned  during  each  time  step.  In  contrast,  DSMC  collision 
procedures  either  give  roughly  linear  scaling  of  computational  expense  with  time  step  size,  or  permit  the  occurrence 
of  nonphysical  collisions  for  which  the  collision  probability  is  greater  than  one.  Alternate  DSMC-based  approaches 
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such  as  collision  limiters21  typically  require  many  times  more  operations  under  these  conditions  than  procedures, 
such  as  those  presented  here,  which  involve  direct  sampling  of  particle  properties  from  known  distributions. 

VI.  Conclusions 

A  method  has  been  presented  for  the  simulation  of  rarefied  gas  flows,  based  on  a  probabilistic  solution  to  the  ES- 
BGK  model  of  the  Boltzmann  equation.  Methods  described  by  Macrossan6  and  Gallis  and  Torczynski8  were 
extended  to  allow  for  the  consideration  of  a  gas  mixture  involving  diatomic  species,  and  a  procedure  has  been 
outlined  to  enforce  exact  momentum  and  energy  conservation. 

Results  from  an  ES-BGK  simulation  were  compared  with  DSMC,  NS  and  experimental  data  for  a  nozzle  flow  of 
the  type  found  in  cold-gas  spacecraft  thrusters.  Generally  good  quantitative  agreement  was  observed  between  ES- 
BGK,  DSMC  and  experimental  data  for  several  flow  properties,  and  the  ES-BGK  flowfield  solution  was  shown  to 
incorporate  thermal  nonequilibrium  effects  which  were  consistent  with  DSMC  results.  Based  on  this  test  case,  we 
find  the  overall  accuracy  of  the  ES-BGK  method  to  be  superior  to  a  traditional  NS  implementation  but  less  than  that 
of  DSMC  for  the  simulation  of  a  rarefied  expansion  flow  involving  a  simple  diatomic  gas.  In  general,  the  proposed 
method  was  shown  to  retain  the  accuracy  of  a  numerical  solution  to  the  NS  equations  in  the  Knudsen  number  range 
where  these  equations  may  be  assumed  valid,  and  at  higher  Knudsen  numbers  accounts  for  effects  of  moderate 
thermal  nonequilibrium  that  are  neglected  in  the  standard  NS  formulation. 

While  the  ES-BGK  method  is  built  on  a  series  of  near-equilibrium  approximations  that  may  produce  large  errors 
in  the  simulation  of  highly  rarefied  flows,  it  avoids  the  calculation  of  individual  molecular  collisions  which 
generally  makes  DSMC  impractical  for  calculations  involving  very  low  Knudsen  number  flows.  As  particle 
movement,  boundary  conditions  and  sampling  procedures  are  the  same  in  the  ES-BGK  method  as  in  DSMC,  the 
method  is  ideally  suited  to  implementation  in  a  hybrid  code  for  the  simulation  of  flows  involving  a  range  of 
Knudsen  number  regimes.  This  avoids  much  of  the  complexity  inherent  in  NS/DSMC  hybrid  codes,  and  reduces  the 
sensitivity  of  computed  flow  characteristics  to  the  interface  location. 

An  optimization  of  simulation  parameters  to  balance  accuracy  and  efficiency  is  out  of  the  scope  of  this  study,  but 
the  lack  of  individual  collision  calculations  in  the  ES-BGK  method  should  allow  for  computational  efficiency 
significantly  greater  than  DSMC  in  near-equilibrium  regions  where  a  relatively  large  time  step  can  be  used.  The 
overall  efficiency  will  however  be  comparable  to  DSMC  in  regions  with  a  higher  degree  of  nonequilibrium,  where 
assumptions  underlying  the  method  are  likely  valid  only  if  the  time  step  size  is  restricted  to  the  local  mean  collision 
time.  This  method  therefore  holds  promise  as  a  compromise  in  both  accuracy  and  efficiency  between  NS  and 
DSMC,  and  may  be  a  possible  alternative  to  NS  solvers  in  a  hybrid  code  with  DSMC,  or  to  continuum  methods 
based  on  higher-order  moments  of  the  Boltzmann  equation  for  the  simulation  of  transition  regime  flows  where 
assumptions  underlying  the  NS  formulation  may  break  down  but  DSMC  is  prohibitively  expensive. 
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Fig.  1  Time  variation  in  temperatures  for  Fig.  2  Grid  geometry  for  ES-BGK  and  DSMC 
rotational  relaxation.  simulations. 
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Fig.  3  Mach  number  contours  for  ES-BGK 
(top)  and  DSMC  (bottom). 


Fig.  5  Density  variation  along  the  nozzle 
centerline. 


y/Rw 


Fig.  7  Density  variation  along  a  radial  plane 
at  x/R,=18.7. 


Fig.  4  Density  contours  for  ES-BGK  (top)  and 
DSMC  (bottom).  Values  are  normalized  by  the 
stagnation  density. 


Fig.  6  Variation  of  rotational  and 
translational  temperature  along  the  centerline. 


Fig.  8  Temperature  variation  along  a  radial 
plane  at  x/Rt=18.7. 
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Fig.  10  Continuum  breakdown  parameters 
along  the  nozzle  centerline. 
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